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ABSTRACT 

We develop an analytical Hamiltonian formalism adapted to the study of the motion of two 
planets in co-orbital resonance. The Hamiltonian, averaged over one of the planetary mean lon- 
gitude, is expanded in power series of eccentricities and inclinations. The model, which is valid 
in the entire co-orbital region, possesses an integrable approximation modeling the planar and 
quasi-circular motions. First, focusing on the fixed points of this approximation, we highlight re- 
lations linking the eigenvectors of the associated linearized differential system and the existence of 
certain remarkable orbits like the elliptic Eulerian Lagrangian configurations, the Anti-Lagrange 
(Giuppone et al., 2010) orbits and some second sort orbits discovered by Poincare. Then, the 
variational equation is studied in the vicinity of any quasi-circular periodic solution. The funda- 
mental frequencies of the trajectory are deduced and possible occurrence of low order resonances 
are discussed. Finally, with the help of the construction of a Birkhoff normal form, we prove that 
the elliptic Lagrangian equilateral configurations and the Anti-Lagrange orbits bifurcate from the 
same fixed point L4. 

Subject headings: Co-orbitals; Resonance; Lagrange; Euler; Planetary problem; Three-body prob- 
lem 
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The co-orbital resonance has been extensively studied for more than one hundred years in the framework 
of the restricted three-body problem (RTBP). In most of the analytical works, the emphasis has been placed 
on the tadpole orbits, trajectories surrounding one of the two Lagrangian triangular equilibrium points, since 
these describe the motion of the Jovian Trojans. Howe ver, the global to pology of the co-orbital resonance 
has been studied in particular by 



Garfinkel 



(1976 19781 and 



Erdi 



(19771, but the interest for the horseshoe 



orbits which encompass the three equilibrium points £3,^4 and L5, remained academic until the discovery 



of the Saturnian satellites Janus and Epimetheus (Smith et al. 1980 Synnott et al. 1981). In Dermott 



and Murray (1981a), general properties of the tadpole and horseshoe orbits are described in the quasi- 



circular case. In particular, asymptotic estimates of the horseshoe orbits lifetime and the relative width 
of this orbits domain are given. But the impossibility to get explicit expressions of the horseshoe orbits 
complicated their study, and theoretical works were replaced by numerical simulations. Thereby, |Christou| 



(2000) showed that the region containing the tadpole orbits is not disconnected from the horseshoe one, 



and that there exist transitions between these two domains. Also, a global study of the phase space of the 



co-orbital resonance was presented in the RTPB by [Nesvorny et al. ( 2002 1 using a numerical averaging of 



the disturbing function over orbital frequencies. Using the same kind of numerical technics, [Giuppone et al.| 



(2010) studied the stability regions and families of periodic orbits of two planets locked in the co-orbital 
resonance. Besides the Lagrangian triangular configurations, where the three bodies in Keplerian motion 



occupy the vertices of an equilateral triangle, these authors found a new family of fixed points (equilibrium 
in the reduced average problem, but quasi-periodic with two fundamental frequencies for the non-average 
problem in inertial reference frame) that they called Anti-Lagrange orbits. Both Lagrange and Anti-Lagrange 
families can be seen as a one-parameter family of stable fixed points, parametrized by the eccentricity. As 
shown in the Fig. 7 of Giuppone et al. (20101), when the eccentricity is equal to zero, the corresponding 



configurations of each family seem to merge in the well known circular Lagrangian equilateral configuration. 
Finally, when the eccentricity increases, the stability regions surrounding these orbits become smaller, the 
one associated to the Lagrangian configuration being the first to vanish. 

In this paper, we will develop a Hamiltonian formalism adapted to the study of the motion of two planets 



in co-orbital resonance. We modify the methods presented in Laskar and Robutel ( 1995 1 in order to get an 



analytical expansion of the planetary Hamiltonian averaged over an orbital period, meaning averaged over one 
of the planetary mean longitudes. This expansion, which is a power series of the eccentricities and inclinations 
whose coefficients depend on the semi-major axes and on the difference of the planetary mean longitudes, 
generalizes the expressions obtained in the RTBP framework by iMorais (1999 20011. Moreover, this one 



containing only even terms in the eccentricities and inclinations (see Sectional), the planar circular motions 
are conserved along the solutions. These quasi-circular co-orbital motions are modeled by an integrable 
Hamiltonian depending only on the semi-major axes and the mean longitudes difference. Contrary to the 



integrable approximations derived by Yoder et al. (1983) or Morals (1999 2001), our model possesses five 
fixed points. Three are unstable and correspond to the Eulerian collinear configurations denoted by Li, L2 
and L3 where the two planets revolve on circles centered at the Sun. The two others correspond to the 
circular Lagrangian equilateral configurations L4 and L5 mentioned above, which are linearly stable if the 
planetary masses are small enough ( Gascheau||1843 ). 

In Section [Sj the linear differential system associated to infinitesimal variations transversal to the plane 
containing the quasi-circular orbits will be studied. Only the directions corresponding to the eccentricities 
will be considered. First, focusing on the fixed points, we will highlight relations linking the eigenvectors 
of the linearized differential system and the existence of certain remarkable orbits like the elliptic Eulerian 
and the Lagrangian configurations, the Anti-Lagrange orbits and some second sort orbits discovered by 



Poincare ( J1892 ). Then, the variational equation will be studied in the vicinity of any quasi-circular periodic 
solution. The fundamental frequencies of the trajectory will be deduced, and possible occurrence of low 
order resonances will be discussed. 

Section[4]will be devoted to the construction of a Birkhoff normal form in the neighborhood of L4. Beside 
the derivation of the fundamental frequencies of any quasi-periodic trajectory lying in this neighborhood, 
we will prove that the elliptic Lagrangian equilateral configurations and the Anti-Lagrange orbits bifurcate 
from the same fixed point L4. Finally, in Section [5] comments and various approaches for future works will 
be presented. 



2. The average Hamiltonian 
2.1. Canonical heliocentric coordinates 



We consider two planets of respective masses mi and TO2 orbiting a central body (Sun, or star) of mass 
niQ dominant with respect to the planetary masses. As only co-orbital planets are considered, no planet is 
permanently farther from the central body than the other, so the heliocentric coordinate system seems to be 



the most adapted to this situation. Following Laskar and Robutel (1995), the Hamiltonian of the three-body 



problem reads 



with 
and 




(1) 



where r^ is the hehocentric position of the planet j, j3j = momj{mo + rrij)^^ and ^j = GiiTio + rrij), Q being 
the gravitational constant. The conjugated variable of r^, denoted by fj, is the barycentric linear momentum 
of the body of index j. In this expression, Hk corresponds to the unperturbed Keplerian motion of the two 
planets, more precisely the motion of a mass /3j around a fixed center of mass too + ruj, while Hp models 
the gravitational perturbations. If we introduce the small parameter s given by 

s = Maxf!^,^), (2) 

Vtoq Too/ 

one can verify that the Keplerian term of the planetary Hamiltonian is of order e and the other one is of 
order £^ which justifies a perturbative approach. 

The choice of these canonical heliocentric coordinates (fj, Tj) may lead to quite surprising results for 
quasi-circular motions. In particular, the famous Lagrangian relative equilibrium, where the three bodies 
occupying the vertices of an equilateral triangle animated with an uniform rotation, is described in terms of 
elliptical elements by ellipses in rapid rotation. More precisely, f j being not collinear to the heliocentric ve- 
locity of the planet j, the Keplerian motion associated to the unperturbed Hamiltonian f ^/(2/3j) — fij/3j/\\rj \ \ 
is not represented by a circle described with constant angular velocity, but by a rapidly precessing ellipse 
whose eccentricity is proportional to the planetary masses. This phenomenon described in the appendix 
(Section [6]) is similar to the question of the definition of elliptical elements for a satellite orbiting an oblate 



body (see Greenberg (1981)). Except this little drawback which occurs only when the considered motion 
is close to the circular Lagrangian configurations, the canonical heliocentric variables are particularly well 
suited to the study of the co-orbital resonances. 

In order to define a canonical coordinate system related to the elliptical elements (a^, Cj, Ij, Xj, Wj, il^) 
(respectively the semi-major axis, the eccentricity, the inclination, the mean longitude, the longitude of the 
pericenter and the longitude of the ascending node of the planet j), we start from Poincare's rectangular 
variables in complex form (Aj, Aj, Xj, —ixj, yj, —iVj) where Aj — Pj^jjijoj, 



^A/ -^^ \/l-ejexp(iti7j), 



V] = V Aj Y V ^ " ^^*^"^ ~ cos/j)exp(il7j). 

This coordinate system has the advantage of being regular when the eccentricities and the inclinations tend 
to zero. It is also convenient to use the non-dimensional quantities Xj = Xjy/2/Kj and Yj = i/j/ y^2Aj which 
are equivalent to Cj exp(i-cuj) and Ij exp(ir2j)/2 for quasi-planar and quasi-circular motions. 

As we only consider the planetary motions in the vicinity of the circular planar problem, the Hamiltonian 
can be expanded in power series of the variables Xj, Yj and their conjugates in the form 

*^^q^-^(Ai,A2)xf^xf xf xf y^^^yf Ff Ff' ] e'(fci^i+fc2A2)^ (4) 




where the integers occurring in these summations satisfy the relation 

Y,{kj + Pj + Qj - P] -q,)^ 0, 



(5) 



known as D'Alembert rule. This relation corresponds to the invariance of the Hamiltonian by rotation, or, 
which is equivalent, to the fact that the angular momentum of the system is an integral of the motion. 
Remark that we will not use this explicit Fourier expansion in this paper, but the D'Alembert rule will play 
an important role. 



According to Poincare (19051, the expression of the angular momentum in Poincare's variables reads 



c-E= 



E 



/ v^S(y,)y^ A,-|:r:,-P-^ ^ 



V 



-V25R(y,)^A,-|:c,P- 



\yj? 



(6) 



/ 



In order to deal with the co-orbital resonance, an appropriate canonical coordinate system iaM 
9 J , Jj , Xj , -ixj , y-j , -iyj ) with 

0i = Ai-A2, 2Ji=Ai-A2, 
6i2 = Ai+A2, 2J2 = Ai+A2. 



(7) 



Inside the 1:1 mean motion resonance, the angular variable 9i varies slowly with respect to 02- Conse- 
quently, the planetary Hamiltonian ([T]) will be averaged over the angle 02- 



2.2. The quasi-circular and planar average problem 

2.2.1. The average problem 

In this paper, we only consider the average Hamiltonian at first order in the planetary masses. More 
precisely, we assume that there exists a canonical transformation which maps the initial Hamiltonian H to 



H{0], Jj,Xj, -ixj,yj,-iyj) = Ho{Jj) + i/i(6'i, J,-, Xj, ~ixj,y-i,-iyj) + 0{e^) 



with 



and 



Ho{Ji, J2) — — 



Pff^l 



Pif^i 



2( Ji + J2V 2( Ji - J2)2 



= Hko 



7j , Jj , Xj , IX j , yj , tyj ) 



Hi{ei, J J , Xj , -ixj , yj , -iyj ) 



1 

2^ 



Hp o (f>{ej , J J , Xj , -iXj , yj , -iy.j)de2 , 



(8) 

(9) 

(10) 



where the map </> satisfy the relation (f j , r j ) — (f>{6j , Jj , Xj , —iXj , yj , —iyj ) • If we denote by {6j , Jj , Xj , ~ixj , yj , 
— iyj) the canonical variables associated to the average problem, we remark that J2 = Ai -I- A2 is a first inte- 



gral of H . It is also easy to prove that the quantities ^ yj y^j 



|y,|V2andEjA, 



\yj\ 



are first integrals too. It is possible to take advantage of these first integrals by reducing the problem by 



Nesvorny et al. 



I 2002 I for the 



^ Other coordinates adapted to the co-orbital resonance have been used by several authors (e.g. 

RTBP and |Giuppone et al.| | [2010[ l for the planetary problem), but these systems, that performed the reduction of the angular 
momentum, are singular when the eccentricities tend to zero. 



means of adapted canonical coordinate system as it is the case with the Jacobi reduction in the spatial 
problem (see Robutel (1995) and|Malige et al^(2002)). The reduction can also be achieved in the planar 
problem leading to two degrees of freedom Hamiltonian system depending on two angles: the difference of 



the mean longitudes and the difference of the longitudes of the perihelion (Giuppone et al. (2010)). These 



reductions introducing some technical issues (addition of a parameter, singularities when the eccentricities 
and inclinations tend to zero), we prefer not to reduce the problem. 



The average Hamiltonian ( |10p depending on the mean longitude only by their difference 9i , the rotational 
invariance of the Hamiltonian given by the relation i5h imposes that H is even in the variables Xj and yj and 
their conjugates. As a consequence, the set xi = X2 = yi = 1/2 = is an invariant manifold by the flow of the 



average Hamiltonian (10). More generally, this property holds for any order of averaging. This implies that 
the part of the average Hamiltonian ( |lOl ) which does not depend on the eccentricities and the inclinations, 
namely Hq{9i, Jj) — H{6i, Jj, 0, 0, 0, 0), is an integrable Hamiltonian. It is worth noting that the one degree 
of freedom Hamiltonian Hq, associated to the circular and planar resonant problem, is a peculiar attribute 
of the 1:1 mean-motion resonance. The next section is devoted to its study. 



2.2.2. The integrable part Hq 

After replacing the vectors Tj and f ^ by their expressions in terms of elliptic elements into the planetary 
Hamiltonian (fl]) , an explicit expression of Hq is obtained by suppressing the terms depending on the variables 
Xj , Xj , yj , y J and the fast angle 62 ■ This leads to the Hamiltonian 



Hn 




(11) 



y/a,ia2 \Jo\ + ^2 ^ 2a\a2 cos B\ 



where the semi-major axis Oj depends on the action J\ and the first integral J2. The constant 2J2 = Ai + A2 
being positive, there exists a strictly positive number a such that 

ATMi +/32VM2 /= 



J2 



(12) 



At this point, it is convenient to define a new couple of conjugate variables (0, J) by translating the action 
J\ as 

«. /77T _ «, /TTT 

\=d. (13) 



^ /3l^/Ml -/32VM2 /= , ^ 
J\ = — — va + J, 



It will also be useful to define the dimensionless (non canonical) action-like variable u by the relation 

J = {/3i+ l32)\/nQau with /ig = QrriQ. (14) 

Now, by a substitution of the relations 

IfJ-o 







(15) 



into the expression (11), the integrable average Hamiltonian Hq can be explicitly expressed in terms of the 



(9, J, a), or {9,u,a) for convenience. Note that the expression (15) allows one to interpret the parameter a 
as a mean value around which the semi-major axes oscillate. 



The figure n] reproduces the phase portrait of the integrable Hamiltonian Hq in coordinates {6, u). It can 
easily be expressed in terms of semi-major axes using the expression ( 15 ) or their first order approximation 
Qj — a ^ [2(— l)^+^a(r7ii + TO2)/toj1 u. The upper plot represents the whole phase diagram for mi — mj = 
10^^ and 1712 — fns = 3 x 10^^ and Q = niQ = a = 1, where the masses mj and ms are close to those of 
Jupiter and Saturn expressed in solar mass. This plot is similar to the well known Hill's diagram (or zero- 



velocity curves) of the non averaged planar circular RTBP (see Szebehely (1967)) although the zero- velocity 



curves are not solution curves of the motion. It is also topologically equivalent to the phase space of the 



average planar circular RTBP when the eccentricity of the test-particle is equal to zero (Nesvorny et al 
2002 Morbidelli|2002 1 . The Hamiltonian system associated to Hq possesses five fixed points that correspond 
to the usual Euler and Lagrangian configurations, and one singular point at u = 9 = which corresponds 
to the collision between the planets. The two stable equilibrium points located at 9 — ±7r/3, u — (see 
the next paragraph for more details) represent the average equilateral configurations that we will denote 
abusively by L4 and L5 by analogy with the RTBP. Each of these points is surrounded by tadpole orbits 
corresponding to periodic deformations of the equilateral triangle. This region is bounded by the separatrix 
S3 that originates at the hyperbolic fixed point L3 at = tt, u « 0, for which the three bodies are aligned 
and the Sun is between the two planets and its separatrix. Outside of this domain, the horseshoe orbits are 
enclosed by the separatrix ^2 that originates at the fixed point L2 {9 = and u < 0). This point, as the 
equilibrium point Li, is associated with an Euler configuration for which the two planets are on the same 
side of the Sun. The last domain, centered at the singularity, is surrounded by the separatrix Si connecting 
the Li point {9 = and u > 0) to itself. Inside this small region, the two planets seem to be subjected to a 
prograde satellite-like motion, the one revolving the other one clockwise. By an enlargement of this region 
( —0.1 < 9 < 0.1), the second plot of Fig. IT] (middle box) shows the splitting of the two separatrices Si 
(red) and ^2 (blue) when the planetary masses are different. On the contrary, for equal planetary masses, 
the phase portrait becomes symmetric with respect to the axis u = 0. It turns out that the equilibrium 
points L^jLijL^ lie on the axis of symmetry, and that the two curves Si and ^2 merge together giving 
rise to a unique separatrix connecting Li to L2- The bottom plot of Fig. [l] describes this phenomenon for 
TOi = 1712 = 5 X 10""*, the other parameters being unchanged. 

A way to estimate the locations of the equilibrium points is to use an asymptotic expansion of the 



Hamiltonian (11) in the neighborhood of u = 0. Two cases have to be considered. The first one arises in a 
domain which excludes a suitable neighborhood of the collision (the distance between the planets has to be 
of order unity) . The second case concerns a small domain enclosing the singularity. In the first situation. 



our goal can be achieved by an expansion of the Hamiltonian (11) in the neighborhood u = 0, assuming 
that the condition 9 = 0(1) is fulfilled. We will see later that this is satisfied in the tadpole region. This 
condition also holds for the horseshoe orbits which do not approach too much the singularity. Denoting by 
r the quantity \/2 — 2cos0, and using the notations ai = mi + 1712, (j'l = mi — TO2 and CT2 = 7ni7n2, the 
expansion of Hq can be written as 

Ho{9, u) = Ga-' [71 + J2{9) + 0{e') + {62(9) + 0{e')) u 

+ {Ti+T2{9)+0{e'^))u^ + u'^R{9,u,e)\ , 

where the coefficients ^k^^k, 5k are given by 

27i = -moai, 272 = ^2(2 - P^ - 2r-i), 

252 = -<Ji<j'i{l-Tf{l + 2T-^), 2ti = -3moafa^\ (17) 

T2 = ala2' [i^l - 3a2)(4 - PV2 - P^^) + 2a2p-3] , 

and the remainder i? is a periodic function of 9 depending on u and of order e. The location of the 



fixed points L3, L4 and L5, as well as the eigenvalues of the associated linearized system, can be easily 



deduced from the expansion ( 16 ). The location of the two elliptic fixed points L4 and L5 is approximated by 
{6,u) — (±7r/3,0 + 0(£^)), which leads to ai = a(l + 0(e^)),a2 = a(l + 0{e'^)). The quadratic expansion 



of the Hamiltonian Hq in the neighborhood of L4 (change 7r/3 in — 7r/3 for L5) is equal to 

2^ 



H, 



"'^^ =-li (f ("^"'^^ - ^'^^ + ^"^)"' + h i'-l 



(18) 



where only the dominating terms in e are retained. Moreover, the frequency associated to this elliptic fixed 
point reads 



i/Q = no 



27 0-1 
4 Too 



1 



0-2 



2?71oCri 



+ 0(£') 



(19) 



where uq — /Iq a ^' ^ plays the role of an averaged mean motion. The location of the hyperbolic point L^ is 



obtained by translating u by a quantity u^^' which cancels the linear term in u in the expansion ( 16 ) when 
9 = 7:. We get the approximation 



,(3) 



(to-i — TO,2)to,iTO,2 

3?no(TO,i + to,2)^ 
which gives in terms of semi-major axes: 

{-ly TOfe mi - m2 



-O(e'), 



Qj = a\ 1 



0{e' 



with j ^ k. 



(20) 



(21) 



3 TOq TTii + ni2 

As a consequence, the quadratic expansion of iJo in the vicinity of £3, whose coordinates are (tt, w--'^')^ reads 

3 5 (al 



H, 



(2) 

0,-^3 



7 
2 a \(T2 6 



,(3))2_^^^(^_^)2 



(22) 



The domain including the tadpole orbits is bounded by the separatrix 53 . The size of this domain can 
be estimated by difi^erent ways. A simple manner to achieve this goal is to calculate the quantity f/3 which 
is the maximal value taken by the action u along 53. By solving the equation Hq{tt/3,U3) = Ho{tt,u^^') 



where Hq is approximated by (16), we get the expression: 

V2a2 



U.= 



\/Srrioaf 



+ 0{e') 



72: 



7TT-1TO2 



•\/3mo(TOi -|-TO2)3 



O(e^). 



(23) 



Similarly, in the perpendicular direction, we can also estimate the quantity Q3 corresponding to the minimal 

value of 6 along ^3 by an approximation of the positive root of the equation Ho{Q^,0) = i/o(7i',M^'^^). The 

solution is given by 

/2 — 1 
63 = 2arcsin(^-- — ) + 0{e) w 23.9°, (24) 



which is a classical result in the case of the RTBP ( GarfinkeT||1977 ) 



The approximation ( 16 ) is not valid for the Euler points Li end L2, the latter being located at a distance 



of order e^'^ of the singularity. In this case, we can use the asymptotic expansion of Hq{0,u), valid since 
u = 0{e°') with < a < 1, given by 



HoiO,u) = ga- 



7-ii^r' 



^2 



I'l 



^'2+Oie^) + {S'2+Oie^))u 



(25) 
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Fig. 1. — Phase portrait of the Hamihonian Hq in coordinates (0, u). The upper box shows the whole space 
for Too = 1, TO-i = mj and TO2 = ms, the parameters G and a being equal to one. The separatrix that 
originates at L3 (53) is plotted in green, while ^2 is the blue curve and iSi the red one. The middle plot is 
an enlargement of the region surrounding the collision point while the bottom plot shows the merging of the 
two separatrices Si and 52 when the planetary masses are equal. Here, their values are toi = 7712 = 5 x 10^^. 



where the coefhcients 7^ , Tk , Sk are 

27-1 = -crr^^27 271 = 271 = -TO-oCTi, 

2S2 — 3aia[, 2ti = 2t[ — — StooCTiCT^ 
4 = 4{al - 3a2)ala^' 



72 =cr2, 



(26) 



and R' is a periodic function of 6 depending on u and of order e. At this accuracy, the two Euler points Li 
and L2 are symmetric with respect to the hne u = 0, and if we denote by u^^' (resp. u^'^') the u-coordinate 
of Li (resp. L2), we have u^^^ = — u^^) + ©(e^/^) and 

Let us mention that the quantities u^^' , u^"^' and u^^' can also be considered as roots of a polynomial equation, 



as it is the case for Euler's configurations in the full three-body problem (see Marchal and Bozis ( 1982 ) or 



Roy (1982)) 



As for the tadpole orbits, the width of the horseshoe region along the u axis can be deduced from the 
equation Hq{tt/3, Ui) — H(){0,u^^^) that is 

C/i = 2-i/26i/6mo-'/Vr'/'a2 + 0(e'/'). (28) 

The minimal angular separation between two planets in horseshoe orbit is solution of the equation Hq{9^^\ 0) = 
i/o(0,u(^^), and is equivalent to 

^ ^ "'^ +0(£2/3). (29) 



3 \6mo 

We note that at this degree of accuracy (neglecting the terms of order e^/^), the two separatrices Si and ^2 
merged in a single curve. 



The equations (231 and (28) allow one to retrieve the result on the relative size of the tadpole and 



horseshoe regions obtained by|Dermott and Murray (1981a) in a very different way. Indeed, we have 



^^2^('!^^i±^V^' = 0(en (30) 

As a consequence, the lower the planetary masses are, the larger the horseshoe region is (with respect to the 
width of the tadpole region) . 

It is worth to mention that, although the average system modeled by the Hamiltonian Hq provides a 
faithful representation of the topology of the problem, it only reflects poorly the dynamics in the domain 
bounded by Si containing the singularity. 

The simplest argument that points out this problem comes from the computation of the orbit frequency 
surrounding the collision point. Indeed, close to the singularity, the Hamiltonian can be approximated by 

— Qmira2ar [aj + )^ ' with a = {mim2lJ.oo)^ . (31) 

It turns out that the frequency of the trajectory that originates a,t 6 = 6q and J = is equivalent to 

y/mim2 no 



Too 



^n^ 



(32) 



As a product of an averaging process, this frequency would be small compared to ng, but it tends to infinity 
when 0Q tends to zero. 

To conclude this section, we will compare the average Hamiltonian Hq to a classical approximation of 
the co-orbital resonance. This model, represented by the Hamiltonian Ha which reads 

3gmo{mi+m2f 2 , Gmim2 ( ^ ^ \ /QQ^ 

Ha = —- u -I cosy = , (33) 

2 a TOim2 a \ V2-2cos0y 
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Fig. 2. — Comparison between the model Hq and its approximation Ha- The level curves of Hq are plotted 
in red while the level curves of Ha are in green. The two approximations fit correctly when the trajectories 
do not come too close to the singularity (left panel) . The right plot shows the inconsistency between the two 
models in the Hill region. 



has been used by Yoder et al. (1983) to study the dynamics of the co-orbital satellites Janus and Epimetheus 



(see also Sicardy and Dubois ([2003i)) . A similar Hamiltonian is also employed to model the 1:1 mean motion 
resonance in the RTBP (see Morals (1999)). This approximation is a particular case of the expression (16) 



obtained by expanding Hq in power series of u and e. For small enough values of u, and if 6 is not too 
close to or to 27r, the Hamiltonian Ha provides a good approximation of Hq in the tadpole region and for 
horseshoe orbits providing u « e^/^. As one can see on the left plot of the figure[2J although the differential 
system associated to the Hamiltonian Ha possesses only three fixed points, its trajectories are very close to 
those of the average Hamiltonian Hq, when they do not approach the collision. This is especially true for 
the tadpole orbits and the moderate amplitude horseshoe orbits. On the contrary, as shown in Fig. [2] (right 
panel), the second model is not valid in the Hill region. But we have to keep in mind that though in this 
region the topology of the average problem is well described by Hq, it is not the case of its dynamics. 



3. Variational equations in the neighborhood of the quasi-circular problem 



3.1. The variational equations 



It has been shown in the previous section that the manifold Xj = yj = is invariant by the flow of the 



average Hamiltonian ( 10 ). In order to study the (linear) stability of this invariant manifold in the transversal 
directions (xj^yj), we have to calculate the variational equations associated to this invariant surface. These 



equations, corresponding to the linearization of the differential system associated to the Hamiltonian (10) in 
the neighborhood of the plane Xj = yj =0, can be derived from the quadratic expansion in eccentricity and 



inclination of the average Hamiltonian H. This expansion can be written in the form Hq + H2 + i?2 with 



(h) _ 



A^) 



i/f ^ = gmim2 (A;,XiXi + B^X^X^ + B,,X,X2 + ^/.^a^a) 



and 



H^"' = Qm^rn^ {A^Y^Y ^ + 5,^1^2 + B^Y^Y^ + A^Y^Y 2) 



(34) 



(35) 
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where the coefficients Ah , Bh , A„ and B„ read 



Ai, 



f aia2 1 

V A3 ^aia2 



COS0, B„ 



V«i«2 A3 



0102 \ iS 

^ 1 



^,, = ^ (aia2(5cos20 - 13) + A{al + al) cos9) - """^^ 



8A5 



2^/aIai"' 



-2i0 



(36) 



B, 



aia2 



aia2 (e 



27^10^ 16A5 
A = W a^ + flg — 2aia2 cos 9 



-3»e _^ gg.e _ 2Q^-^e^ ^ g(^2 ^ aj)e-^'') 



The formulas (34), (35 1 and (36) generafize the expansion given by Morais (1999 2001) in the case of the 
elhptic RTBP. 



The variational equations in the vicinity of a solution lying in the plane Xj ^ yj = and satisfying 



IdHn 



, w), ii 



IdHn 



de 



.u) with c =(/?!+ /32)y7loa, 



take the form 



and 



X2 
Y2 






iQmim2 



ArM, 



A2 B^ 



A2 ^-u 



Xi 
X2 

Yi 
Y2 



Mh{0,u) 



MviO,u) 



Xi 
X2 



Yi 
Y2 



(37) 



(38) 



(39) 



where 9 and u are deduced from the solutions of the equations ( 37 1 , and the Aj implicitly depend on u 



by the relations (15). As these solutions are periodic (except if their initial conditions are chosen on the 
separatrices Si to ^3) the linear equations (38) and (|39[) are periodically time-dependent. As a consequence, 



their solutions cannot generally be expressed in a close form. A notable exception occurs at the equilibrium 
points of the system (37). Indeed, here, the variational equations become autonomous and consequently 



integrable. Then we will first begin to study these special cases. Before going further, let us mention that in 



this paper, we will not study the "vertical" variational equation (39). Indeed, due to its strong degeneracy. 



the study of this linear equation is not sufficient to understand the local dynamics. To this aim, the use 
of higher order terms of the Hamiltonian is necessary (at least the forth degree in yj). To be convinced, 
it is enough to look at the matrix Afi,, defined in (39), for 9 ~ 7r/3 and u — 0. Indeed, at £4 the matrix 



vanishes and the quadratic Hamiltonian does not provide any information about the dynamics in the yj 
directions. Then, this situation requires a careful analysis of the structure of the Hamiltonian in order to 



deal with potential bifurcations in the vertical direction, as it is pointed out by Jorba ( 2000 ) in the case of 
the bicircular problem. Therefore, we postpone this study to a future work. 



3.1.1. Dynamics around the fixed points 



For the equilateral configurations {9 — ±7r/3), neglecting the quadratic terms in e, the matrix Mh takes 
the following expression 

,27 no f 1712 —77126* 



Mh 



mo 



-TTiie 



TOl 



(40) 
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This matrix possesses two eigendirections associated to the eigenvectors 



Vi=[ '''""' ) and V2=[\ ), (41) 



-mi 



whose eigenvalues are respectively 



. 27 mi + ma 

Vi — —I no and V2 — U. (42j 

8 mo 

These eigenvectors have a precise physical meaning. Along the neutral direction, the one which is collinear 
with V2, the two eccentricities are the same and the angle Aw = tui — W2 separating the two apsidal lines 
is equal to 7r/3 at L4 and — 7r/3 at L5. These configurations clearly correspond to the Lagrangian elliptic 
equilibria, which are fixed points of the average problem, and consequently of the linearized average problem 
at L4 or L5. This is the reason why the associated eigenvalue V2 vanishes. Along the direction Vi, the orbits 
satisfy the relations 

ai = 02 = a, 6 = ±7r/3, miei = 772262, and Atu = tui — 'uj2 = 6 + 11. (43) 

This corresponds to an infinitesimal version of the Anti-Lagrange orbits found numerically by |Giuppone 



et al. (20101. On these trajectories the elliptic elements ai, 02, ei, 62 and 9 are constant. Only the two angles 



zui and W2 precess with the same frequency equal to 

27mi+m,2 , , 

ffi = wi = -z no, (44) 

8 m-o 

in such a way that the angle Aw is constant. As a consequence, this family of periodic orbits is transformed, 
after reduction by the rotations, in a family of fixed points, which is exactly what have found [Giuppone et al.| 



(2010) in the reduced problem. Of course, the family that we found along the eigenvector Vi of the linearized 
system provides only an infinitesimal approximation of the Anti-Lagrange family in the neighborhood of L4, 
but we will show in Section [4] that this linear approximation can be generalized to any degree using Birkhoff 
normal form. 



3.1.2. The Eulerian fixed point L3 

By evaluating the matrix Mh{9,u) at {9,u) = {ir^u''^^) and neglecting the terms in e^ and more, the 
matrix of the linearized system at L3 reads 

K = *^— f""' "0- (45) 

8 mo V ™i "^1 / 

This matrix possesses two eigendirections associated to the eigenvectors 

"■' - ( L": ) "^ ''i - ( -1 ) ■ (*) 

whose eigenvalues are respectively 

, .7mi+TO2 , 

Wi — I- tiQ and V2 = 0. (47) 

8 Too 
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As in the equilateral case, the direction V2 corresponds to the unstable Euler configurations where the two 
planets are in the two sides of the Sun (the eccentricities are equal and the perihelia are in opposition). 
The other direction is more interesting. In this case, the perihelia are in conjunction and the eccentricities 
verify the relation mici — m2e2. As for L4, the method developed in Section 4.1 makes possible to prove 
the existence of a one-parameter family of periodic orbits that bifurcates from L3 and is tangent to V( at 
this point. Remark that, at least close to L3, this family has been numerically computed by [Hadjidemetriou 



et al. (2009), and were previously found by Poincare (18921 as a solution of second sort |^ (see 



Chenciner 



3.1.3. Euler L I and L2 equilibria 

As the Hamiltonian Hq does not reflect properly the dynamics in the neighborhood of the collision 
between the two planets, the linearized problem at Li or L2 will not be considered in the present paper (see 



the end of Section 2.2.2) 



3.2. The general solution of the variational equation 

Now, let us study the general case. This corresponds to writing the variational equation ( [SS] ) around 
a periodic solution of frequency i^. According to the Floquet theorem (see Meyer and Hall (1992)), the 
solutions of the variational equation take the form 



z{t) = P{iyt)eTq){At), 



(48) 



where A is a constant matrix and P(V') is a matrix whose coefficients are 27r-periodic functions of tp- As, if 
Z is a fundamental matrix solution to the variational equation along a 27r/z^-periodic solution, one has the 
relation 



Z{t + 2TTiy-^) = Z{t) exp {2t:v-^A) 



(49) 



and the solutions stability of the variational equation depends on the eigenvalues of the monodromy matrix 
exp (27rz/^^ A) . As a consequence, if we start the integration at i = from the identity matrix, after a period, 
we get the relation exp {2'kv~^A) = Z{2'kv~^). Thus, this matrix and its eigenvalues can be deduced from a 
simple numerical integration of the variational equation. If the eigenvalues modulus of this matrix are equal 



to one, the solutions of ( 38 1 are quasi-periodic. Moreover, their fundamental frequencies are v, 51 , 52 , where gi 



and 32 are equal to the eigenvalues arguments of the monodromy matrix multiplied by v /{2'k). Fig. p^ shows 
the results corresponding to planetary masses mi = mj and 7712 = rns. From the numerical computations 
of the monodromy matrix eigenvalues, we conclude that solutions of the variational equation are always 
quasi-periodic, and thus the invariant manifold of the quasi-circular orbits xi = 2:2 = is transversally 
stable, at least in the directions associated to the eccentricities. This property seems to hold for every value 
of planetary masses that we have tested, that is nii — TO2 = 10~^ and m2 = 0.3 x nii = 0.3 x 10^^ with p 
ranging for 3 to 8. Figure [3] shows the behavior of the frequencies ly, gi and 172 along a section of the space 
phase. The red curve corresponds to v, the green one to gi and the blue one to .92- The initial conditions 



are chosen on the segment 6 — tt/3 and < m < u^^\ u^^'' being defined by the relation (27 1 as the positive 



intersection of the line 6 = n/S with the separatrix Si. This plot shows clearly two different dynamical 



^Methodes nouvelles de la inecaniquc celeste Vol I, Chap III, &47: "Solutions de la seconde sorte" 
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domains: the inner one filled with tadpole orbits ranging from m = to u^^', and the outer domain for 
u*^^) < u < M^^-* populated by horseshoe orbits. 

Inside the inner region, the libration frequency v decreases from the value v^ ~ -^27(7711 + ?7i2)/(4mo)no 
~ 0.0936 yr~^ to zero when the separatrix ^3 is reached. The frequency gi associated to the precession of the 
periastra evolves smoothly between 27{mi+ni2)no/{8mo) ~ 0.00439 yr~^ at L4 and 7(toi + m2)no / {8mo) ~ 
0.00114 2/r~^ approaching 53. The box located in the upper left corner of the plot details its evolution for 
< M < 0.0055. As shown in this figure, the last frequency 32 is always very small with respect to the 
other ones. It starts from zero and reaches zero again at the separatrix, being at least always twenty times 
smaller than gi. Because in the tadpole region, the frequency v is of order ^/e and gi of order e, these two 
frequencies do not generate low order resonances (some of these resonances are indicated by vertical black 
dotted lines), except in a very narrow neighborhood of 53. As i' tends to zero at ^3, in both sides of this 
separatrix, the two curves intersect and v becomes smaller than gi. Using the estimates derived by|Garfinkel| 



( 1976 1978 ) in the RTBP, one can easily show that the frequency v reaches a logarithmic singularity where 
it tends to zero as — (log|M — u^'^^l)"^. Consequently, the slope of the curve associated to v is very steep 
and then the low order resonances occur only very close to the separatrix, in a region which is intrinsically 
unstable. 



0.14 



0.12 



0.08 - 



0.06 



0.02 - 



0.0055 

0.005 - 
0.0045 

0.004 
0.0035 

0.003 



0.001 0.002 0.003 0.004 0.005 0.006 



0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 



Fig. 3. — Variations of the fundamental frequencies of the variational equation, u is plotted along the X-axis 
and the frequencies, in rad/yr, along the Y-axis. The red curve represents the evolution of the libration 
frequency v along the segment 9 = 7r/3, < u < u'^', while the green (resp. blue) curve is associated 
to 5i (resp. (72)- The vertical black dotted lines indicate the location on some resonances between gi and 
ly. The box plotted in the upper left corner of the figure is an enlargement showing the behavior of gi for 
< M < 0.0055. 



The situation is more interesting inside the horseshoe domain. Indeed, if 52 remains always very small 
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with respect to the frequency gi , this one increases significantly after the crossing of the separatrix ^3 . Note 
that this behavior is aheady mentioned by'Morais|(|1999) where the figure l.a of her paper dedicated to the 



RTBP corresponds approximately to our figure [3] for u ranging from to about 0.008. After this value, gi 
keeps to increase, reaching the resonance v — 2gi at u sa Q.Ol and even the 1:1 resonance between v and gi 
approaching the end of the horseshoe domain materialized by the separatrix Si . After the crossing of the 
two curves representing the frequencies v and gi, gi remains temporarily above v. The situation is reversed 
quickly as v tends to zero when Si is reached. 

The behavior of the three fundamental frequencies v, 51 , 52, described for mi = mj and TO2 = ms, seems 
to be very weakly mass-dependent. Indeed, the simulations performed with the mass sample presented above 
converge to the same conclusions. First, the frequency 32 is always small with respect to gi and of course to 
p. Second, no significant resonance which may destabilize the average system occurs in the tadpole region, 
excepted in a narrow area surrounding S3: the lower the planetary masses are, the larger the ratio I'/gi is. 
And third, in the horseshoe domain, low order resonances involving 1/ and gi always occur, in particular, 
the 1:1 close to Si is crossed two times. If these low order resonances generate chaotic behaviors in the 
average problem, it is not necessarily this mechanism that dominates in the full (non averaged) three-body 
problem for planetary masses comparable to those of Jupiter or Saturn. Indeed, [Laughlin and Chambers] 



(2002) deduced from numerical simulations of the planetary three-body problem that horseshoe orbits, even 
starting with the two planets in circular motion, are unstable for planetary masses satisfying the empirical 
relation (mi -I- m2)/{mQ + mi + 7712) > 0.0004. This limit corresponding approximately to two Saturn's 



mass planets around the Sun. One can find comparable simulations in different cases in Dvorak (2006). In 



a nice paper by Barrares and Olle (2006), this behavior is studied more carefully. In the case of the RTBP, 
these authors prove that the invariant manifolds associated to L3 deeply penetrate the region populated 
by horseshoe orbits, generating a large chaotic region, whose size increases with the mass of the secondary 
(the mass of the primary being fixed) . They even mentioned the possible heteroclinic intersections with the 
invariant manifolds associated to Lyapunov orbits around Li and L2- A similar mechanism probably acts 
in the planar planetary problem. As suggested by C. Simo (private communication), not only L3 invariant 
manifolds, but also invariant objects of periodic orbits existing in the vicinity of the previous manifolds are 
supposed to be involved in the process, as it is the case in the spatial RTBP. This phenomenon, acting in short 
time-scale, plays a major role in the instability of the horseshoe regions, even for zero initial eccentricities. 
For moderate to small planetary masses, the portion of the horseshoe orbit region intersecting the invariant 
manifolds mentioned above shrinks to a narrow region, excluding transitions between the L^ region and the 
neighborhoods of Li and i2- In absence of short time-scale chaos, the destabilizing effect of the resonances 
involving the frequencies v and gj can dominate, at least locally, the dynamics of the full problem, as it is 
the case for the average one. 



4. Beyond the quadratic approximation: BirkhofF's normal form and family of periodic 

orbits 

Let us begin with the study of the dynamics in the neighborhood of L4 (the discussion would be the same 
at L5). We first start with the linearized system at this point, or equivalently, with the quadratic expansion 



of the average Hamiltonian in the vicinity of L4 . Using the notations (17) and ( 34 1 , this expansion takes the 
form 



''\2^„...^^fTW 



mio-^r+v2u'+m"\ (50) 
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the coefficients 771 and 772 being deduced from (18). A symplectic diagonalization of the associated Hamil- 
tonian system allows one to define a new canonical coordinate system (zq, zq, zi, zi, Z2, ^2) that reduces the 

previous quadratic form to 

2 

K2^^1qZqZq. (51) 

As for the variables Xj and Xj, the coordinates Zj and Zj are linked with the relation Zj = —izj. L4 being 
an elliptic equilibrium point, the coefficients 7^ are purely imaginary. More precisely, we have 70 = iv, 
71 = igi and 72 = ig2 = 0. Consequently, if < j. A:, / < 2 are three distinct integers, the set defined by the 
equation Zj = z^ := is a one-parameter family of periodic orbits of the linearized system parametrized by 
the complex number zj. The frequency, which is given by |7i|, is the same for every orbit of the family. Let 
us denote J^q the family parametrized by zq that corresponds to the quasi-circular motions (ei —62 — 0). 
The one parametrized by zi corresponding to the linear approximation of the Anti-Lagrange orbits will be 
denoted J"{. And the last one, governed by Z2, which contains the Lagrangian elliptic configurations, will be 
symbolized by 7^2- 

Let us now consider the term of degree greater than two in the expansion of the average Hamiltonian 
in the neighborhood of L4, and let us write this expansion as 

K = K2 +Y,Kp with Kp= Y, -i^zl°~zfzfifzf~zf, (52) 

p>3 qeX'e.p 



where 



n-l 

2?2„,p = {q = (50,90, • • • ,'Zn-i,gn-i) e N^"/ |q| = ^ {\qj\ + \qj\)=p}. (53) 

J=0 



All the coefficients q are not allowed in the summations (52): as \xi\'^ -\- |x2p — |zip -I- |z2p the D'Alembert 
rule is still valid in coordinates Zj , and the non-zero coefficients 7q verify the relation qi + (72 = 9i + 92 ■ This 
last relation imposes that the total degree of the monomials z\^zf'z'^z'^ is even, thus, the manifold given by 



the equation zi = Z2 = is still invariant by the fiow of the Hamiltonian K defined in ( 52 1 . It turns out that 



the family J-q is not only an invariant set of the linear problem ( 51 ) but also of the full average Hamiltonian 



(52). This statement also holds for the family F2 including Lagrange's configurations. Indeed, as we know 
that these configurations exist as fixed points of the average problem and that we always have 6 = 7r/3 and 
fli = 0,2 (or u = 0), Zq = along the family J-2. In addition, the relations An7 = 7r/3 and ei = 62 = constant 



impose that zi — 0, according to Section 3.1.1 This implies additional constraints on the coefficients of the 



Hamiltonian K. Indeed, as every element of this family is an equilibrium point, the Hamiltonian K fulfills 
dK _ 9K 

dzj dzj 



the conditions f^ = ^ = when |zo| = |zi| =0. This is equivalent to the cancellation of the coefficients 



of the terms 

{z2~Z2r, Z^{Z2~Z2Y\ ~Z^{Z2~Z2Y\ ^1^2(^2^2)«^ ~Z1Z2{Z2~Z2Y^ . (54) 

As regards the Anti-Lagrange family, the relations |zo| = |z2| =0, which characterizes its infinitesimal 
approximation T\ does not hold. Indeed, if these previous relations are preserved by the linear flow of the 



system associated to ^2, it is no more the case by the fiow oi K. If the Lyapunov center theorem (see Meyer 



and Hall (1992)) could be applied to K, it would show the existence of a one parameter family of periodic 
orbits originating at L4 and tangent to J-"{, whose periods would be close to 27r/|7i| in the neighborhood of 
L4. Unfortunately, the coefficient 72 being equal to zero, the hypothesis of the latter are not fulfilled. To 
overcome this difficulty, we use a more elaborated method, based on the construction of a Birkhoff normal 
form. 
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As mentioned above, the use of the coordinates (zj,Zj), provides an elementary parametrization of 
the famihes J^q and J-2. It is possible to build a coordinate system (CjiCj) for which the Anti-Lagrange 
family possesses the same kind of parametrization than the two other families, that is |Col = IC2I — ^^nd 
C,i depending on the element of the family. This coordinate system can be chosen among one of those that 
reduce the Hamiltonian K to its Birkhoff 's normal form. In this context, the Birkhoff transformation consists 
in the construction of canonical transformations that act on homogeneous polynomials of given degrees in 
order to eliminate non-resonant monomials. These are the monomials which are not of the form 

^qo ~qo qi ~qi q2 ~q2 /'rr\ 

Zq Zq Z]^ z-^ Z2 Z2 . yoo) 

More precisely, this transformation is performed iteratively, each step being dedicated to the normalization 
of a given degree. An elementary transformation is defined by the time-one map of the flow of an auxiliary 
Hamiltonian w„ defined as a solution of the equation 



where $„ contains non resonant monomials of degree n (see Morbidelli (20021)). This equation being linear. 



it can be solved monomial by monomial. The resolution of the equation ( 56 1 for ^„ = Zq" Zq° z\^ zf^ z'^ z^ 



introduces the divisor 70 (go — go) + 7i(gi — gi)- If we assume that 70 and 71 are rationally independent, 
which is generically the casa^ the denominator cancels only if go — go and gi = gi independently of qi and 
qi. Using the D'Alembert rule, the only monomials involving divisors equal to zero are z'^ z^^ z^ z^ z'^ z^ 
which are not eliminated from the Hamiltonian. Consequently, the Birkhoff normal form can be computed 
at any degree. Let us denoted by (CjiCj) the normalizing coordinates. By construction, the coordinates 
(zj,Zj) and (CjXj) ^^^^ related by expressions of the form: Zj — Q + 02{CjXj) with Q — ~iCj- Then the 
Hamiltonian reduced to a Birkhoff normal form reads 

2 

N{QXi) = Y.^qCpCp+ E 7;(CoCo)'°(ClCl)'^C2C2)''^ (57) 

9=0 qo+qi+q2>2 

where the 7' are complex numbers such that the coefficients of the monomials (C2C2)^^ vanish. As an 
example, the Birkhoff normal form corresponding to mi — mj and TO2 = ms computed up to the fourth 
degree in Q , Q reads 

-0.093622 <oCo - 0.00439 iCiCi - 2450.55 Co Co + 472.218 CoCoCiCi 
+253.10 C0C0C2C2 - 38.0734 Ci'Ci - 1.17035 C1C1C2C2, 

where only a few digits of the coefficients are given here. Remark that the "linear" fundamental frequencies, 
namely the coefficients of the monomials iCoCo and jCiCii are negative real numbers. As it is more convenient 
to deal with positive quantities, we have decided to change the sign of these frequencies in the previous 
sections. In the coordinates (Cj,Cj), the Hamiltonian system associated to N is trivially integrable. In 
particular, its phase space is foliated in 3-dimensional invariant tori carrying linear flows. In other words, 
using the angle-action variables {(pj , Ij ) defined by the relations Q — \/lj^^'^^ 1 one can verify that the actions 
Ij are integrals of the motion, and that every solution is quasi-periodic with fundamental frequencies equal 



^As 71/70 ~ %/27{mi +m2)/mo/4:, only high order resonances can occur. This allows one to build the normal form up to 
a high degree, typically of order 1/y/e. In our numerical application, the first potential small denominator involves terms of 
degree 48. 



- 18- 

to ujj = ^. Among these solutions, we focus now on those that are members of the famihes Tj, that is the 
solutions satisfying the relations Ik — Ii — {j,k,l pairwise distinct), or equivalently Cfe = = 0- Using 
the transformation that reduces the Hamiltonian K to its normal form up to a given degree, say 2n, and 
taking into account the symmetries of the transformatiorFl one can show that the families are parametrized 
as follows. The family J^q containing the quasi-circular periodic orbits is given by 

^0 - Co + /(Co, Co), zi^z2^ 0, Co - -*Co e C, (59) 

/(Co, Co) being a polynomial of degree 2n in (Co, Co) whose lower order terms are quadratic. The family Ti 
associated to the Anti-Lagrange orbits reads 

^0 = /'(CiCi), ^1 = Ci + CiQ(CiCi), Z2 = Ci^(CiCi), (60) 

where P, Q and R are polynomials of a single complex variable of degree n whose lower order term is of degree 
one. Of course, we still have Ci = — ^Ci G C. As mentioned above, the elliptic equilateral configurations J^2 
are still given by: 

zo = zi = 0, Z2 = C2, C2 = -«C2 £ c: . (61) 

J^i is the most interesting of these three families. Indeed, the quasi-circular family J^) is well known and 
its orbits are already represented in the figures [T] and [2J The elliptic equilateral configurations J^2 has also 
been extensively studied since their discovery by Lagrange. In addition, its expression in terms of elliptic 
elements is well known since it corresponds to ai — a2, ei ^ 62 and 9 = Azu = 7r/3 (or — 7r/3 for the family 



starting from L5). On the contrary, the family Ti has been only partially studied by Giuppone et al. (2010) 



and ,Hadjidemetriou and Voyatzis (2011) 



The use of a Birkhoff normal form allows one to get any desired information concerning this family, 
providing that the orbits of the family are contained inside the domain of validity of the normal form. 
Practically, this is not the case for the whole family, at least some portion of J^i including L4 is contained 
in such a domain. In order to estimate this region, we have computed the relative difference between K and 
N along J^i using the expression 

. > N ^ \K{zo, zi,Z2,zo, zi, Z2) - N{0, Ci, 0, 0, Ci, 0)| 
|iV(0,Ci,0,0,Ci,0)| 



where the values of Zj are deduced from Ci by the relations (60 1. As p(0) = 0, and in order to be consistent 
with the approximations done during the computation of the average Hamiltonian, we consider that the 
normal form is relevant while p(Ci) < e^- We have estimated that the thirtieth degree was a good compromise 
between the precision of the normal form and its number of terms. Once defined this domain in which the 
normal form is relevant, a linear transformation allows one to express the Zj (deduced from Q) in terms of 
6,u,xi,X2 and to deduce the expression of J-"i with the help of the elliptic elements. This is shown on the 
left panel of the figure |4] in the particular cases rrii — mj, 7712 = mg. Fig. |4ja displays the evolution of the 
eccentricity ei versus 62 along the family J^i (red curve). The maximal value of 62 for which the condition 
p(Ci) < £^ (here p(Ci) < 10~^) is fulfilled is 62 = 0.23, which corresponds to ei « 0.066. Let us note that, 
we have p(Ci) < 3 x 10^^^ as long as 62 < 0.12 and that the precision obtained using the Birkhoff normal 
form is comparable to the machine epsilon. In this domain, ei seems to depend linearly on 62, the slope 
of the (red) line being equal to "012/ rni. The difference between the green curve, which shows the variation 



It is not necessary to detail this transformation, but the key point lies on the fact that it takes the form C,j 

f 



fj{zi, Z2, Z3, zi, Z2, Z'i) where the polynomial fj possesses the same symmetries as 
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Fig. 4. — Evolution of the elliptic elements and of the fundamental frequencies in function of e-i , along the 
families J-\ and J-^- Left panel: elliptic elements. e\ red curve, panel (a). Q in red and Atu (green) in panel 
(b). Panel (c): ai — 1 (red), a^ — \ (green). Right panel: the fundamental frequencies i/, g\ and gi are 
represented in (d), (e) and (f). The green curves correspond to frequencies computed along J^i, T2 in red. 



of TOiei — 771262 versus 62, and the dashed black line (ei = 0), indicates that the relation miCi 
fulfilled only at the origin of the family T\. 



771262 IS 

Hadjidemetriou and Voyatzisl (|2011|) suggest that along this 



family, e\ and 62 tend simultaneously to one regardless of the planetary masses. Fig. [4]b shows how the 
angles d in red and Atu — 180° in green move away from their value at the origin when 62 increases. Basing 



on numerical simulations, Hadjidemetriou and Voyatzis (20111 suggest that the angles Q and Atu tend to 



180°, when the eccentricities tend to one, which would correspond to a triple collision. The last figure of 
the left panel. Fig. |4]c, shows the slight deviation of the semi-major axes from the equality oi = 02 = 1. 
Practically, ai — 1 is plotted in red, while the green curve corresponds to a^ — 1. This figure shows that, 
at least for 62 < 0.23, the variations of the semi-major axes are very small (of order e^) compared to the 
other elliptic elements. The situation may be different for large values of the eccentricities, but this is not 
mentioned in the literature. 

Remark that, with the help of the analytical expression of J^i, we analyze a relatively small portion of 



the family T\ compared to the region studied numerically in Giuppone et al. (2010) and Hadjidemetriou 



and Voyatzis (2011) where the eccentricities reach 0.8. In contrast, our analytical study allows us to access 
to more information. First, it provides a complete understanding of the dynamics of all quasi-periodic 
trajectories lying in the validity domain of the Birkhoff normal form. Second, using an analytical expansion 
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of the eigenvectors of the differential system (38), we can estabhsh rigorously that, at the beginning of the 
family T\, the orbits satisfy the relation m\e\ — 771262, which has been empirically deduced from numerical 



simulations in Giuppone et al. (20101. Third, it allows us to compute straightforwardly the fundamental 



frequencies associated to each trajectory belonging to a given family. Indeed, for Ti^ the derivative of the 
normal form iV with respect to /; is the frequency of the corresponding periodic orbit of the family (this 
frequency is zero in the particular case of J-2). The normal frequencies are obtained by derivation with 
respect to the two other action variables. These three frequencies are plotted in Fig|4jd-f for the families T\ 
and J^2- The fundamental frequencies associated to the family Tq are not represented here for the simple 
reason that the normal form furnishes the same values as in figure [3] at least in a neighborhood of the circular 
equilateral configuration L4. The frequency v (resp. 51, g-i) is plotted in Fig. [ij-d (reps, [ij-e, |4]-f). The red 
curves correspond the equilateral family Ti while the green curves are associated to T\. 

Although these frequencies and their derivatives are equal at the origin of the families, their behaviors 
along T\ and J-2 are very different. As shown Fig. [4]-f, the frequency gi is obviously equal to zero all 
along the Lagrange family since these trajectories are fixed points of the average problem. On the contrary, 
computed on the family J-\ this frequency increases to a (local) maximum although it remains small in the 
considered interval. According Fig. |4]-e, g\ changes only very slightly for the equilateral family, but very 
much for T\. Remark that the quantity 27: / gi, which seems to increase with the distance to L4, is the period 
of the orbits belonging to J-i. Regarding u (Fig. |4]-d), the frequency associated to J-i seems to reach a local 
maximum, while the one corresponding to T2 increases. 

What can be said concerning the behavior of the fundamental frequencies outside of the validity domain 
of the normal form? One thing is clear about the equilateral configurations: when their eccentricity increases, 
a critical value depending on the mass ratio (momi +17101712 + 71211712)/ {mo + mi + 1712)'^ is reached, leading to 
a period-doubling bifurcation where the family looses its stability ( Roberts]|2002 Nauenbe"rg||2002 1 . Conse- 
quently, for 62 > 0.23, the frequency v is supposed to keep increasing, until it reaches the resonance 2v = n, 
where n is the planetary mean motion (close to one ii a = Q — mo = 1). This is certainly the mechanism 
that was acting when Giuppone et al. ( 2010| ) observed the shrinking of the stable region surrounding the 
equilateral equilibrium, and finally its fading when the eccentricity grows. The way that the family J^i ends 
is less clear. In fact, at high eccentricities, only numerical simulations of these orbits have been performed 
( Giuppone et al.|2010 Hadjidemetriou and Voyatzis|2011 ), and when e does not exceed 0.8. Hadjidemetriou 
and Voyatzis (20111 suggest that for high eccentric orbits, the two eccentricities coincide, and that and 



Atj7 tend to tt. This would imply that the Anti-Lagrange family Ti, and the Euler family originating at L3 
intersect, or end at a triple collision. This conjecture has to be checked. 



5. Concluding remarks 

In this paper, we developed a Hamiltonian formalism adapted to study the motion of two planets in co- 
orbital resonance. This analytical formalism intends to unify several works dedicated to the 1:1 mean-motion 
resonances like the formulations developed by Erdi (1977) or 



Morals 



( 1999^,2001) in the case of the RTBP, 



but also models obtained by Dermott and Murray ( |1981a ) and Yoder et al. (1983) aiming to understand the 
dynamics of the two Saturn's satellites Janus and Epimetheus. 

Our approach consists on an expansion of the average Hamiltonian in power series of both planetary 
eccentricities and inclinations. To make the study of the tadpole orbits as well as the horseshoe orbits 
possible, an expression of the mutual distance valid for all values of = Ai — A2 has been introduced in 
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the Hamiltonian. Contrary to the other authors who modeled the distance between the two planets by the 
term \/2 — 2 cos 0, wc have chosen to introduce the divisor \J a\ + a^ — 2aia2 cos 9. This changes drastically 
the topology of the integrablc problem associated to ei — £2 = Ii ^ I2 = 0. Indeed, the usual model, 
which possesses three fixed points corresponding to £3, L4 and L5, is singular when 6 — 0^ regardless of the 
planetary semi-major axes values. Our approximation gives rise to two additional fixed points corresponding 
to the Euler points Li and L2- The singularity, that is usually identified to a line in the usual model, 
is here reduced to a single point that corresponds to the collision of the two planets in the same circular 
orbit, that is ai = 02,^ = 0. Thus, the topology of the two problems is very different. Indeed, with the 
first approximation, the phase space is divided in three distinct regions: two symmetrical libration regions 
around L4 and L5 respectively, and a third one, populated with horseshoe orbits that encompass the three 
equilibrium. Inside this last region, the semi-major axes tend to infinity when the angle 9 approaches zero, 
which is obviously not very realistic. With the average model presented in this paper, the two regions 
surrounding L4 and L5 are practically the same as in the usual model, while the horseshoe region bounded 
in a domain lying between the separatrix emanating from L3 and the one originated from L2. This model 
can be useful to simulate captures or transitions between different kinds of trajectories under the influence 
of weak dissipations, or slow migrations. Indeed, contrary to the usual model, the non-resonant region is 
better separated than the resonant horseshoe region. 

For small eccentricities, the global topology of the problem is similar to the one described in |Nesvorny| 



et al. ( 2002J ) in the RTBP framework, using numerical averaging methods which are not limited to moderate 
eccentricities and inclinations. Although we are constrained by the size of eccentricities and inclinations, 
our model possesses at least two advantages. On the one hand, this average problem, as long as the number 
of terms of its Hamiltonian is not too large, allows fast numerical simulations using large time-steps. On 
the other hand, the present analytical formulation of the problem can help to obtain theoretical results 
concerning the stability inside the co-orbital resonance. If much has been done in the vicinity of the equilateral 



equilibrium points, especially in the RTBP (see Gabern et al. (20051 and references therein), the theoretical 



stability of horseshoe orbits remains an open problem. 

With the help of a Birkhoff normal form, we have shown how the equilateral family T2 and the Anti- 
Lagrange family J^i bifurcate from the circular equilateral configuration L4. If the behavior of the family 



J-2 is well known from its beginning at L4 to its termination by a period-doubling bifurcation (Roberts 



2002), the same cannot be said for the family J^i. At this point, we only have conjectures concerning the 



termination of this family. This might be a triple collision, and could be related to the end of the Euler aligned 
configurations originated at L3. A similar question, which is not discussed in the present paper, concerns 



the so-called quasi-satellites family (see Hadjidemetriou et al. (2009); Giuppone et al. (2010)) which could 
also end by collisions when the eccentricities tend to one (an alternation of two kinds of double collisions 
involving on the one hand. Sun and a first planet, and on the other hand, the second planet and the Sun). 



A last point should be mentioned. In Section |3.1[ the vertical variational equation has been set aside 
because the quadratic part of Hamiltonian in inclination was equal to zero. A careful study of this situation 
would reveal interesting bifurcation phenomena giving rise to families of remarkable orbits, as in the case 



of the RTBP (Perdios and Zagouras 1991 Marchal 2009) or in the general three-body problem with equal 
masses (Chenciner and Fejoz||2008). Finally, a lot remains to be done in that field. 
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6. Appendix: L4 in heliocentric canonical elliptic elements 

Let us assume that the three bodies describe a circular Lagrangian equilateral configuration where p is 
the length of the triangle sides. The heliocentric coordinate system can be chosen such as Vj — pUj where 



cos ipj 

sin ipj 




with ifi — Lot and Lp2 ~ ojt -\- 



(63) 



the angular velocity cj of the relative equilibrium satisfying the third Kepler law uj"^ p^ = fi = Q{mQ+mi+'m2). 
The elliptic elements {aj, Cj, Vj, Wj) can be derived from the canonical heliocentric coordinates (fj, r^) using 
the relations 

2/3, llrjil 2a J ' 



K, 



E, 






COS tUi 



(64) 
(65) 



and 



e/Ej-Uj- 



-^ = 7 rj rfc , with 7 = 1- 

a straightforward computation leads to the expressions 



^, {j,k)e{l,2} and jV^, 



2p 



and 



1% 



mj. 



Pj 

1 A^ /3fe 



A^ 



mo m5 



I ) 7-^ 



TTiQ + ?Tij 2 /ij TTiQ y "* Pj rno \ 2 rriQ 



1 Ufc. 



(66) 
(67) 

(68) 
(69) 



According to (68), the semi-major axis of the planet j is a time-independent quantity approximated by the 
expression 



aj= p[l 



nik mi + 7712 
7710 rno 



+ Oie^) 



(70) 



which is slightly larger than the radius p of the configuration. As Ui • U2 = 1/2, the expression (69) shows 
that the eccentricity (modulus of Ej) is constant, and that the ellipse rotates with an angular velocity equal 



to w. A first order expansion of (69) gives 



E, 



and 



777o V 2 / 

\/3mk 



2 7n,o 



+ Oie^). 



We deduced from ( 66 ) that the true anomalies Vj of the planets satisfy 



cos Vi 



2v377io 



+ 0{e^). 



(71) 
(72) 

(73) 
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